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We describe a numerical code that solves Einstein's equations for a Schwarzschild black hole in 
spherical symmetry, using a hyperbolic formulation introduced by Choquet-Bruhat and York. This is 
the first time this formulation has been used to evolve a numerical spacetime containing a black hole. 
We excise the hole from the computational grid in order to avoid the central singularity. We describe 
in detail a causal differencing method that should allow one to stably evolve a hyperbolic system 
of equations in three spatial dimensions with an arbitrary shift vector, to second-order accuracy in 
both space and time. We demonstrate the success of this method in the spherically symmetric case. 
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I. INTRODUCTION 

A key goal of numerical relativity is to determine the gravitational radiation produced by the inspiral and coalescence 
of two black holes in a decaying binary orbit. The importance of solving this problem is heightened by the possibility 
that gravitational wave detectors such as LIGO, VIRGO, and GEO may observe gravitational waveforms from binary 
black hole coalescence within the next decade. Not only could a comparison of measured waveforms with numerical 
simulations provide a crucial strong-field test of general relativity, but in addition, accurate templates produced by 
these simulations could significantly increase the sensitivity of waveform measurements and reduce the uncertainties 
in astrophysical parameters derived from gravitational wave data. 

Although constructing a numerical simulation of binary black hole coalescence is a difficult and unsolved problem, 
several recent advances have brought us closer to a solution. One key advance is the development of so-called 
apparent horizon boundary condition (AHBC) methods which treat black holes by evolving only the regions 

of spacetime that lie outside apparent horizons. These methods take advantage of the fact that information cannot 
emerge from within the apparent horizon of a black hole (which, assuming cosmic censorship, is contained within the 
event horizon). Without AHBC schemes, the spacetime singularity that inevitably forms inside a black hole eventually 
causes numerical simulations to terminate, typically on a time scale of order 10 to lOOM, where M is the mass of the 
system. For this reason, AHBC methods may be crucial for solving the binary black hole problem, where one hopes 
to evolve two holes long enough to see them orbit, coalesce, and eventually settle into a final state containing a single 
(Kerr) black hole. 

Another promising development is the construction of manifestly hyperbolic formulations of Einstein's equations 
|p[-|T3| . Not only do these formulations offer insights into the mathematical structure of general relativity, but they may 
also be better-suited for numerical solution than the usual ADM equations, which are not manifestly hyperbolic. 
The reason for this is twofold: First, there exists an extensive literature concerning stable and efficient numerical 
methods for solving hyperbolic systems of equations jl^] . Some of these methods have been successfully applied to 
hyperbolic formulations of general relativity by Bona and Masso Q . The second reason is that a non- hyperbolic set 
of equations can present a fundamental difficulty for black hole simulations that employ an AHBC approach. 

To understand this difficulty, consider a non-hyperbolic set of equations, or even a hyperbolic set that has charac- 
teristics lying outside the local light cone. Suppose such a system is to be solved on a domain that includes a black 
hole. Although the physics guarantees that nothing can emerge from the hole, the equations do not know this, and 
nonphysical (gauge) modes can propagate outward through the apparent horizon. Solving such a system of equations 
on a restricted domain that excludes the interior of the hole is mathematically well-posed only if appropriate boundary 
conditions are imposed on the horizon. While it should be possible to impose explicit horizon boundary conditions 
to fix the coordinate system [^ , it is unclear which boundary conditions are appropriate for dynamical variables. 
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particularly in the general three-dimensional case in which one may have a nonspherical horizon and a significant 
amount of gravitational radiation. 

Now consider a hyperbolic set of equations with characteristics that never lie outside the local light cone. In this 
case, future-pointing characteristics inside the apparent horizon (which must be an outgoing non-timelikc surface) can 
never intersect the horizon itself, so that quantities at or outside the horizon cannot depend on the interior region. 
Consequently, one can solve this set of equations on a restricted domain that excludes the interior of the hole, and 
one can do so without imposing boundary conditions on the horizon. 

In this paper, we concentrate on the hyperbolic formulation of Einstein's equations originally proposed by Choquet- 
Bruhat and York , hereafter referred to as the CBY system. This formulation has several advantageous features. 
First, the characteristics of this set of equations are extremely simple: they lie either along the light cone or along the 
normal to the current time slice. This guarantees that no information, not even gauge information, can propagate 
acausally. Second, the equations admit an arbitrary shift vector, and the characteristics are independent of the choice 
of shift. Finally, the fundamental variables in this formulation are spatially covariant three-dimensional tensors. These 
tensors directly measure spacetime curvature, and from them one can form all components of the spacetime Weyl 
tensor. 

While we believe that the CBY formalism holds considerable promise for numerical simulations, particularly when 
combined with AHBC methods, there is no experience with solving this particular set of equations on a computer. 
Therefore, before expending the significant effort required to implement the CBY equations in a full three-dimensional 
code, it is important to demonstrate that such an approach is feasible. 

Accordingly, we have developed a numerical code that solves the spherically symmetric Einstein equations using 
the CBY formalism. We evolve a Schwarzschild spacetime using a causal differencing scheme, and we use an AHBC 
method to avoid the central singularity. This is the first time the CBY equations have been used to evolve a numerical 
spacetime containing a black hole. 

Our code runs in parallel, and is rigorously second order convergent. As described in detail elsewhere [ pT| , it can 
run for times in excess of lOOOM provided certain constraints are regularly enforced. Our code is based on the DAGH 
1^ software package originally developed for the Binary Black Hole Grand Challenge Alliance. It is written in C-I--I-, 
and uses Fortran-90 numerical kernels. The DAGH system contains support for adaptive mesh refinement, but we 
have not yet taken advantage of this feature. 

Our code provides an important demonstration that the CBY formulation works well in numerical simulations. It 
allows us to study the details of implementing the CBY equations in a simple setting, and it provides a testing ground 
for apparent horizon excising schemes and causal differencing algorithms. It also serves as an important check on a 
code that solves the CBY equations in three spatial dimensions, a code that is currently under development and will 
be described elsewhere. 

We employ results and methods specific to spherical symmetry as little as possible, so that our techniques are readily 
generalizable to the three-dimensional case with Cartesian coordinates. For example, we do not use logarithmic radial 
coordinates [|[^] despite the great advantages they provide for spherically symmetric codes. Likewise, although there 
are many shift conditions that work well with AHBC methods in spherical symmetry js) , we consider only those that 
can be applied in the general three-dimensional case. 

Recent progress has also been made in spherical symmetry by choosing a global coordinate system that has desirable 
properties near the horizon [Q. The spatial gauge used by depends on the concept of an areal radius, and is thus 
applicable only in spherical symmetry. We forego such a coordinate choice in favor of a more general approach. 

In section^ we summarize the CBY formulation of Einstein's equations, and we specialize this formulation to the 
case of spherical symmetry. In section III , we present a key ingredient of our code: a causal differencing scheme that is 
second-order convergent and has a stability criterion that is independent of the shift vector. We describe this scheme 
in detail for the general case of three spatial dimensions, and then apply it to the spherically symmetric case. In 
section |^we describe the AHBC method employed at the inner boundary of our grid, and the shift conditions that we 
use in order to implement this method. These shift conditions not only allow us to control the motion of the apparent 
horizon through the grid, but they also prevent coordinate singularities that may result from differential stretching 
or compression of the grid in the remainder of the spacetime. In section ^ we discuss the boundary condition that 
we impose at the outer boundary of our grid. This boundary is ideally at spatial infinity, but in practice it is placed 
at a la rge b ut finite radius. In section Vl we present the results of rigorous convergence tests using our code. In 



section 



VII, we close with a short discussion of our results. 



II. EQUATIONS 
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A. The CBY Formalism 



Here we summarize the fundamental equations and variables used in the CBY representation of general relativity. 
For details of the CBY formulation and a derivation of the equations, see Q. 
We write the metric in the usual 3+1 form 



(1) 



where N is the lapse function, /?' is the shift vector, and gij is the three-metric on a spatial hypersurface of constant 
t. 

Define the variables 



A (In TV), 
Djtti. 



(2a) 

(2b) 
(2c) 
(2d) 
(2e) 
(2f) 



Here D is the three-dimensional covariant derivative compatible with the three- metric gij , the time derivative operator 
is 



and £ denotes a Lie derivative. The quantity Kij is the usual extrinsic curvature. 
If we assume that the time coordinate satisfies the harmonic slicing condition 



Ut = 0, 



then the vacuum evolution equations take the form 

9nfl„ = -2NKu 



doN 



(3) 



(4) 

(5a) 
(5b) 
(5c) 
(5d) 
(5e) 



di^Lij — 
doMkij 



-NX, 



N 



doau - NDjUoi = N 



(5f) 

(5g) 
(5h) 

(5i) 



where H is the trace of ify , ^(y) 



{Xij + Xji)/2, the quantities F'^y are the connection coefficients associated with 



the spatial covariant derivative D^, and 

J,j = g,j [H{L -H^ + a'^ak + a^) 

+ K'''{4HKm - 2Lm - AKmkKjr - 2aka, - 2afc,)] 
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- H{3L,j + Gaiflj + iaij + lOKikKj) 

+ afe(4Af(y/ - SM^j) - 4a(jM,)fc'= (5j) 
'h^ + a'aj + 2a] - 2K^^K.jk) . (5k) 



+ ai 

Eq. ( pd| ) is the harmonic shcing condition (^. 

There are considerably more variables and equations in the CBY formalism than in the usual ADM formalism. 
However, the form of the equations is much simpler in the CBY case. While the right-hand sides of Eqs. (H) contain 
many terms, these terms consist solely of algebraic combinations of the dynamical variables and involve no derivatives. 
Eqs. (|5|-||) are tensor wave equations whose characteristics are along the light cones. Eqs. (|5a|-]5^ are even simpler — 
they drag information normal to the surfaces of constant t, that is, along zero-velocity (with respect to the normal) 
characteristics. There are no other characteristics in the system. 

Although we have eliminated some gauge freedom in the choice of lapse function by imposing the harmonic time 
slicing condition (|^), the shift is unspecified and completely arbitrary. The shift is not a dynamical variable in this 
formalism, in the sense that it obeys no evolution equation, and that it appears in Eqs. (^ only through the time 
derivative operator d^. Instead, the shift is an auxiliary gauge variable that may be freely chosen on each time slice, 
and may even change discontinuously from one slice to the next. 

In vacuum, the constraint equations include 

Q^L\ + K'^K.j + a'a^ + a\ (6a) 
= M^,, - M,/' (6b) 



= -I- lia.i + Mj/ (6c) 



= R^j - U.J + HK,j - 2K,kK^ - a,aj - a^j, (6d) 



where Rij is the three-dimensional Ricci tensor. Eqs. ( |6a| ) and (6b) are the familiar Hamiltonian and momentum 
constraints rewritten in terms of the CBY variables, and Eq. ( |6c[ ) is a result of harmonic time slicing. Eq. (33) is the 
familiar ADM evolution equation for K^j^ which in the CBY picture becomes a constraint on L^j. Eqs. (2c|), (p^), 
(H), and the usual relation between T'^ij and derivatives of gij are also constraints that must be satisfied at all times. 
All constraints are preserved by the evolution equations. 



B. Spherical Symmetry 



We write the spherically symmetric three-metric in the general form 

(3)ds2 ^ ^2 ^^2 _^ B^r^{d0^ + sin^ 61^02), (7) 

where (r, 9, (f>) are the usual spherical coordinates. Spherical symmetry reduces the number of dynamical variables 
(including the connection coefficients) from 67 to 16. Define 

TrT = 2BrV\r - 2Brr%, = -—T^ee = - „ . ^ ^ (8a) 

Br Br sm 

ar = a\ = a"*^ (8b) 

Lt = L\ = (8c) 

Kt = K\ - K^^ (8d) 

MrT = Mr\ = A//^ (8e) 

Mrr = M\r = M^^r. (8f) 



where the subscript T denotes "transverse" . Then the vacuum evolution equations (0) take the form 
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where 



daA = -NAK\ 
doBr = -NBtKt 

doKrr = NLrr 

SoKt = N{Lt + 2Kt^) 
doN = -N^{K\ + 2Kt) 
doar — Naor 



doUT 



N 



{2MTr - MrT " UrKr 



-aor 



A2 

2Ktclt 



2A^Br 

<9or^T = -N [KrVrT + 2Br{arKT + M^t 
doMrr = N 



doT^rr = \KrrO,r 



^K\{2MTr - MrT " a^-KT) 



2Br \ 



— Lt 



~ N [ — r^rrOOr + ar{arK^r + Mj-rr/A'^ + ttor)] 



^ N 



Qr 



A^ 



r T 

Br 



doarr - N—ao, 
or 

d 

doMrrr ~ N — Lrr = N [(fl^ - 2T''rr)Lrr + 2K'' riKrrUr + Mrrr)] 

d 

doMrT - N — Lt = N \2KT(arKT + 2MrT) + a-riT] 

or 



dfjLr 



N d 
A^d^^ 
N d 



OoLt 7^-^ MrT — 

A-^ Or 



N 
N 



_ J„, _ ^^'rrMrrr ^ £rI(^_/^2 _ ^^Tr) 

A'^ Br 

r' rrMrT ^ rT / , , 

-Jt 'Z_LL _|_ _Ll_{^MrT + M' 



A2 



A^Br 



Tr) 



Jt = (Lrr + arr){KT ~ K%)^ - LtK' 



ar^K'' 



-2aT{K\ + Kt) + 2Kt^ - 2K\Kt'^ 

+2K\^Kt - + ^{AMTr - 3M^t), 
A'^ 

Lrr{hK\ - AKt) + ar^{K\ - IOKt) + 2arr{K\ - 3Kt) 



+Krr{5K\^ - 6K\Kt + 2Kt^) - 



8MrT , 



Qr 



^^{2Kt - K\) + 2MrTK\ + ar Ka^^ + 3a„ - Ur)^ 



+ AuT - 2K\{K\ - 3Kt) 



The constraints m) become 



2Lt + % + 2Kt^ 



2 1 
K\ + 2aT + -rrfflr 
A^ 







aor + ar{2KT + K'' r) + 



MrT - MTr = 



A2 



2MrT = 



(9a) 
(9b) 
(9c) 
(9d) 
(9e) 
(9f) 

(9g) 

(9h) 
(9i) 

(9j) 
(9k) 
(91) 

(9m) 
(9n) 
(9o) 

(9p) 



(9q) 
(9r) 

(9s) 



(10a) 
(10b) 
(10c) 
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1 

1 



2A^Bt 



or 



or 2Br 



Krr{2KT ~ K%) - Lrr = 



1 



KtK\ -gt-Lt^O. 



(lOd) 
(lOe) 



The additional constraints (|^, (|^), @), and the usual relation between T'^ij and derivatives of gij take the form 



dr 



Kr 







M' 



d_ 

dr 

TrT 



Tr 



Kt - MrT = 



2Br 

|-(ln7V)-a 

ay 



- rT 



or 

d 
ar 







■flr = 



(11a) 
(lib) 
(11c) 

(lid) 
(He) 

(111) 
(llg) 
(Uh) 



III. CAUSAL EVOLUTION METHOD 

Here we present the causal differencing method we use to evolve Eqs. (||) from one spacelike hypersurface to 
the next. Straightforward differencing schemes typically become unstable for large shifts, which are needed for 
the implementation of AHBC methods. Our method is second-order accurate and has a stability criterion that is 
independent of the shift vector. We emphasize that our method is not specific to Eqs. (|^), but can be used to handle 
advective terms in any system of first-order evolution equations. Similar causal differencing methods have been used 
by [HH] and [|9| in order to treat large shifts in the standard ADM formulation of Einstein's equations. 

Our causal differencing method is independent of the actual form of the shift vector. We only place two restrictions 
on the shift: First, it must be a smooth functional of the dynamical variables and of the space and time coordinates. 
Second, it can be computed to second-order accuracy, given second-order values for these variables and coordinates. 
The particular prescription that we use for computing the shift will be discussed in section IV. In this section we 
assume only that such a prescription exists. 

Although the code described in this paper assumes spherical symmetry, our causal evolution method is general. 
In this section we first describe our method for the general case of three spatial dimensions plus time, and then we 
specialize to spherical symmetry. 



A. Overview of the Method 

Figure ^ shows an initial spatial hypersurface labeled by < = and a subsequent spatial hypersurface labeled by 
t = to + At. We wish to evolve quantities defined at the point (^0,2:') to the point (<o + At,a;'), that is, along the 
vector t shown in Figure |^. Here t and are the coordinates defined in Eq. (|l|) , and the vector t is given by 

r= Nn + (12) 

where n is the unit normal to the hypersurface t = tg. 

A large shift vector P tends to cause stability problems in most numerical schemes. Some schemes, including many 
implicit ones, are unconditionally unstable whenever t is non-timelike, as is the case in Figure ^ Other schemes can 
be made stable for an arbitrary shift, but only at the expense of a very small time step At. 
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FIG. 1. Spacetime diagram illustrating the relation t — t + P and showing both the {t, a;*) and the {t, S') coordinate systems. 
The vector t must always lie within the light cone. This is not true for the vector t for a sufficiently large shift (3. 



In order to construct a differencing scheme that works for an arbitrary shift, we introduce an auxihary coordinate 
system. First define a new timehke vector 



i=Nn = t- p. (13) 



Then define new coordinates such that 



i = t (14) 
x' = x\x\t) (15) 
£ix'' = 0, (16) 

and such that the spatial coordinates i* coincide with at t = to . The new coordinates and their relationship to the 
vectors t and i are shown in Figure |l|. Partial derivatives with respect to the new coordinates (f, i*) are given by 

d d . d 

dt^dt'^'d^' ^^^^ 

A = ^A. (18) 

dx^ dx^ dxi 

Our method works by breaking up each time step into two substeps, as illustrated in Figure ^: First, we evolve 
quantities along the vector t, that is, we evolve using the {t,x^) coordinate system, from the points on the slice 
t = to in Figure ^ to the points on the slice t = to + At that are labeled by circles. We then complete the timestep 
by interpolating from the p oint s labe led by circles to those labeled by dots. These two substeps will be considered 



separately in sections III C and III D below 



B. Transforming into the {t,x^) system 
Each evolution equation in the system (0) can be written in the form 
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t^to + At 



B 



t = tn 



FIG. 2. Two-dimensional spacetime diagrams illustrating the two coordinate systems used in our causal differencing scheme. 
Solid and dotted arrows denote the normal vector t and the time vector t at each grid point on the initial time slice. Solid dots 
represent grid points with particular values of a;', and circles represent grid points with particular values of The dots and 
circles coincide at t = to, but not at t = to + At. This is because is constant along t while is constant along t. Diagrams 
A and B represent the same spacetime. The only difference is that A is drawn in the computational coordinate system (t,^*), 
where the time axis lies along t, and B is drawn in the normal coordinate system {t, i'), where the time axis lies along t and is 
normal to the spatial slices. 



doT - Q^S' = R, (19) 

where T, the quantity being evolved, is not necessarily a scalar, but may be a coordinate-dependent object. We wish 
to rewrite Eq. (|l9[), which is defined in the computational coordinate system (t, a;*), in terms of the normal coordinate 
system (f, J*). 

If we consider all quantities to be defined in the (i, x*) basis, we can rewrite the da operator in the {t, i') coordinate 
system as follows: 

do = £{ 
= £t — £ p 

- ^ 
~ dt 

d d d 

= (20) 
at ax^ 



so that Eq. (n9[) becomes 



where 



£ 







(22) 
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(23) 



Here we have used Eqs. (^3|), and (|l7|). The first two hnes of Eq. (|2^) are coordinate independent, but in the 
third line we have assumed the {t,x^) basis in order to write £t — d/dt. In the fourth hne we have separated the 
Lie derivative along /3 into two pieces: The advective piece, P^d/dx^ , is the one responsible for the instability that 
often arises when one tries to evolve along t with a large shift vector. This piece is eventually absorbed into the time 
derivative d/dt. The remaining piece, when operating on some quantity T, describes the change in T induced by 
the change in basis vectors along [3. The operator Cp vanishes when operating on a scalar. Furthermore, CpT does 
not actually contain any derivatives of T, but only contains derivatives of (3. Therefore, no spatial derivatives of T 
appear in Eq. (pl|). 

Note that Eq. ( |2l| ) is not the same as Eq. ( p^ ) transformed into the (t, 5*) basis. By splitting the Lie derivative 
along (3 into two pieces, we have derived an equation that describes the evolution of T defined with respect to the 
(t,a;*) basis along a path of constant i*. The coordinate system used to define tensor components, (i,a;'), is different 
from the coordinate system used to label spacetime points during the evolution, (t, i'). 

Note also that we have introduced an additional auxiliary variable: the Jacobian dx^ /dx^ that appears in Equa- 

and (|l8|), we can find the rate of change of the Jacobian along the vector i: 

d_ 

dt 

To compute the Jacobian, we evolve this equation along with Eq. (|2^). Because we set = at the beginning of 
each time step t = to, the initial value of dx^ /dx"^ on each time step is the Kroeneker delta 51- 

Finally, we require derivatives of /3, which appear in the operator and also on the right hand side of Eq. ([24|). 
Assuming that we have some prescription for choosing the shift given the values of x\ t, and the dynamical variables 
at each grid point, we simply use this prescription to compute (3, and then we obtain its derivatives either analytically 
(if the shift is analytic) or by finite differencing. 



tion im). From Eqs. (17 



dx^ 



dx^ dx 8(3 
dx'^ dx^ dx^ 



(24) 



C. Step 1: Evolve along t 

The first step of our causal differencing method is to evolve Eq. ( ^l| ) and the auxiliary equation ( p^ along the 

vector t, from t — t^ to t — to + Ai. Although this can be done using any standard differencing scheme, the algorithm 
described in this section assumes a scheme with two time levels; a three-level scheme (such as leapfrog) requires a 
slight modification of the algorithm. We use a Macormack predictor-corrector method, which we illustrate with a 
simple wave equation in spherical symmetry, written in first-order form in the (t, f) coordinate system: 

~-^Q = R (25) 
dt or 

^Q^-^P = S. (26) 
dt or 

Here R and S are arbitrary functions of P, Q, r, and t. Given a discrete set of uniformly spaced grid points fj, 
we denote P and Q at grid point i and time i = t„ by P" and Q" respectively. To compute P and Q at time 
t = tn+i = tn + At, we first compute initial guesses P and Q using the "predictor" step 

pn+l ^ pn ^ |i (g»^^ - Q^') + R(P^\ Q^) (27) 

or' - + 1^(^+1 - m + s{pp,Q?), (28) 

where Af — ri+i — fi. The quantities P"^^ and Q"^^ are first order accurate in both space and time. Notice that 
the finite difference approximation to the spatial derivative is one-sided. 

Once we have the predicted quantities, we then compute P and Q at time t — t„+i ~ tn -\- At to second order in 
both space and time using the "corrector" step 
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QT^ = \ [q'I^^ + Ql + %iP':^^ - PP-i) + siPr\ Qr+')) • (30) 

The spatial derivatives are taken in the opposite direction to the predictor step. This ensures that the first order error 
terms introduced by the one-sided derivatives in the corrector step cancel those produced by the one-sided derivatives 
in the predictor step. This cancellation would be spoiled by substituting for Q"^^ in the spatial derivative term 

of Eq. (^9|) or by substituting P"+i for _P"+i in the spatial derivative term of Eq. (^0[). However, it is irrelevant for 
accuracy or stability whether the right hand sides R and 5* in the corrector step are computed using predicted values 



of P and Q as in Eq. (|29|) or using corrected values of P and predicted values of Q as in Eq. (30). We use corrected 
values in the right-hand sides whenever possible because it minimizes the memory required for storing temporary 
variables on the computer. 

The above Macormack scheme is stable (disregarding the boundaries, which will be discussed in section ^) whenever 
At < Ar. For a wave equation with characteristic speed v, the stability condition for the above scheme is the familiar 
Courant condition vAt < Af. 

To obtain a corrected value for a particular variable, we require predicted values for all quantities appearing in the 
equation for that variable. For our prototypical equation (|^), in order to compute a corrected value for T, we must 
have predicted values not only for T, Q, S"*, and R, but also for the Jacobian dx^ /dx^ and for /3 and its derivatives 
(which appear in the operator Cp). The predicted value of the Jacobian is obtained by evolving Eq. ( p^ to first order 
using the Macormack predictor step. The predicted value of the shift is obtained by using our shift prescription to 
compute P from the predicted values of the dynamical variables. 



D. Step 2: Interpolation 



One entire numerical time step should take variables defined on discrete points with particular values of x*, and 
compute quantities at the same values of but at a later time. In Figure ||, this corresponds to taking values defined 
at the discrete points on the slice t — tn and computing quantities at the solid dots on the slice t = tn + At. However, 



when we evolve along t as described in section [H C, we compute quantities at the points that correspond to the circles 
on the slice t = tn + At in Figure ^. 

Our causal differencing method therefore requires a second step, namely, interpolation of our dynamical variables 
from the (t, x*) coordinate system back into the {t, x') coordinate system, that is, from the circled points to the dotted 
points on the slice t = tn + At in Figure ||. The values of x^ at the dotted points and the values of i* at the circled 
points are known; both are the same as the values of at the appropriate grid points on the initial slice. To perform 
the interpolation, we must also know either the values of at the circled points or the values of at the dotted 
ones. 

From Eq. ([l7|), the change in the coordinate along the vector i is given by 



(31) 



Therefore, if we evolve Eq. ( |3l| ) along with Eq. (^ll) using the Macormack scheme, we obtain the value of x* at each 
of the circled grid points at t — tn + At in Figure^. This allows us to interpolate quantities from the circles to the 
dots, working in the coordinate system. 

Note, however, that the circled grid points are, in general, not uniformly spaced in x*, as shown in Figure ^A. 
Instead, the dotted grid points are uniformly spaced in x^ While this poses no difficulty for the spherically symmetric 
case discussed in this paper, in the general three-dimensional case interpolating from an arbitrary set of points onto 
a uniform grid is a nontrivial numerical problem that cannot be treated very efficiently. Much easier and much less 
costly in terms of computer time is to interpolate from a uniform grid to an arbitrary set of points. 

To handle this difficulty, notice that if we evolve along the vector t instead of along the vector t, the coordinates x* 
remain constant and the coordinates x* vary. The change in x' along the vector t is given by 



9x* 9x* 



(32) 



where we have used Eq. (|17|). Therefore, if we evolve Eq. (32) along the vector t using the Macormack scheme, we 
obtain the value of x* at each of the dotted grid points at t ^ tn + At in Figure 0. This allows us to interpolate 
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quantities from the circled points to the dotted ones, working in the coordinate system. The circled points are 
uniformly spaced in i', as shown in Figure ^B. We thus interpolate from a uniform grid (in i') to an arbitrary set of 
points. This can be done relatively easily in three spatial dimensions. 

There is a subtlety in evolving Eq. ( |32| ) along the vector t using the Macormack scheme: In order to implement 
the corrector step of Eq. (^), we require predicted values of the quantities [i^ dx^ /dx^ at the dotted points, but these 
quantities are only known at the circled points. However, from the predictor step of Eq. (|3^), we already have predicted 
values of at the dotted points. Therefore, we use these values to interpolate the predicted values of ^^dx^ jdx^ from 
the circled points to the dotted ones. In this case, we are again interpolating from a uniform grid to an arbitrary set 
of points. 



E. Implementation in Spherical Symmetry 

In spherical symmetry, we solve equations (H) on a numerical grid that extends from some value r = Train just 
inside the apparent horizon to a large radius r — Vmax ■ We denote our numerical grid points by , where i runs from 
imin to imaxi Corresponding to the innermost and outermost grid points. We use a zone-centered grid, so that the 
innermost and outermost grid points do not correspond to rmin and rmax- Instead, we set 



where 



' max ' mtn 



''max ^min 



- 1 Ar 



When evolving along i, the spherically symmetric vacuum evolution equations (W take the form 



dt 

dt 

at 
d 

In. 

dt 
d 

d 



dt 



dt 



dt 



Tr 



dr d 

^NAK\ + A——p'' 
or or 

'NBrKT 

dr d 
Or or 

N{Lt + 2Kt^) 
'N'^{K\ + 2Kt) 
Naor 
N 

TrTao 



df d 
dr dr 

{2MTr ~ MrT - arKr)^ 



2A^Br 



2KTaT 



N 



drd_ drd_fdrd_ 

"T i rr r^-P "To Q~ I o 

Or or or or \ or or 



^TrT = -N [KrTrT + 2Br{arKT + M^t)] 
df d 



N 



dr dr 

2MTr{KT + K\) + [MrT + arKT){KT - K\) 
TrT ( Lrr \1 df d 

^2B-r\^-^'^) ^^^^"-d-rdf^ 



(33) 
(34) 

(35a) 
(35b) 
(35c) 
(35d) 
(35e) 
(35f) 

(35g) 

(35h) 

(35i) 

(35j) 
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at or or 



df d J, 
Or Or 



d 

Ot 



-.Mr 



N df d 
dr df 



-rrr^^drr — N 



Qr 

+ aor 



'P'^Gf'p 

— — 

dr df 



Br 



dt 



-.Mr' 



dt 



dt 



^rdf d ^ 

Nj-jrzLrr 

dr dr 



^df d ^ 
N——Lt 
dr dr 



N_dfd_^ 

JSp. Qj. Qf 



df d 
+ 3M,„ — — /?'■ 
dr dr 

N [2KT{arKT + 2MrT) + arLr] 
df d 



dr dr 



d ^ N df d 
dt^-^d-rdf^"-^ 



N 



N 



l2 + "b7 



Mr 



2MTr 



df d J, 
dr dr 

-Jt- 



rT 



(MrT+MTr) 



(35k) 

(351) 
(35m) 
(35n) 

(35o) 
(35p) 



where the right hand sides (pq|-ps|) are unchanged, and we have included the Cp terms explicitly. Note the second 
derivative of (3^ in the equation for T^rr- This term results from applying £^ to a non-tensorial quantity. Because 
/?'" is not an unknown in the system of equations (^5|) , but is instead an auxiliary variable that may be, for example, 
given analytically as a function of the coordinates, this second derivative should not spoil the hyperbolicity of the 
system. 

5l]), and (p2) become 



In spherical symmetry, Equations 



df df 



dt dr 
dt 



d_ 

dt 



df 
dr 



df\^ d(3^ 



dr I df 



(36a) 
(36b) 

(36c) 



Eqs. ( p6| ) provide values for df /dr to be used in the corrector step, as well as coordinate information for the interpo- 
lation step. 

We use cubic interpolation for causal differencing. This is accurate to fourth order in Af. Quadratic interpolation 
would also suffice, but linear interpolation would not yield a scheme that was second order convergent in time. The 
reason is that linear interpolation makes an ©(Af)^ error on each time step, so that after TV time steps the error 
is of order Af. This is because for a fixed total time itotai that we wish to evolve, the Courant condition requires 

7V-itotal(Af)-l. 



IV. SHIFT VECTOR AND INNER BOUNDARY CONDITION 

In this section we describe how we choose a shift vector /3, and how this choice affects how we handle the inner 
boundary of our computational domain. 



Because both the CBY formalism and the causal differencing scheme discussed in section |I| place no restrictions 
on we are free to choose any shift we wish. Although setting (3 — Q \& the simplest choice, it is often useful to 
employ a nonzero shift vector. One technique is to use the shift to simplify the form of the Einstein equations, for 
example, to eliminate particular components of the spatial metric pO| ] . The disadvantage of this approach is that it 
involves an actual change in the equations being solved. Instead of Eqs. (|^), one would be solving a different (but 
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physically equivalent) system of equations that would include the shift as a dynamical variable, and would no longer 
be hyperbolic. 

We instead use the shift for a different purpose: to allow us to truncate our computational domain just inside 
the apparent horizon, so that we evolve only the exterior region. We thus avoid the spacetime singularity inside the 
black hole. If we were to attempt such a truncation with (3 = 0, numerical grid points originally located just outside 
the apparent horizon would soon fall in, and any grid points located inside the apparent horizon would eventually 
encounter the singularity. This is because for zero shift, numerical grid points follow the world lines of normal 
observers, and these world lines are necessarily timelike. 



A. Shift at the Inner Boundary 

We force the inner boundary of our grid, r = Tmim to hover within a grid spacing of the apparent horizon. To 
accomplish this for a static black hole, we choose (3 near the horizon to point along the outward normal to the 
horizon, and we choose its magnitude so that the local coordinate speed of light in that direction is zero. The horizon 
is then approximately stationary with respect to the spatial coordinates. For the spherically symmetric case, the local 
coordinate speed of light in the outgoing radial direction is given by 

dr N 

so we set at the horizon equal to N /A. 

We track the apparent horizon at each time step, and retain only the grid points that lie on the outside. This is 
done via a masking algorithm that labels grid points outside the apparent horizon as valid, and those inside as invalid. 
Invalid points are never used in the computation, just as if those points did not exist. Because we use a cell-centered 
grid, the inner boundary r = rmin is located half a grid spacing inside the innermost valid grid point. 

Numerical errors typically cause the horizon to drift through the grid, even if one tries to lock it in place by forcing 
/3'" = N/A. While this drift seems to cause minimal difficulty with either the stability or accuracy of the code, for 
coarse-resolution runs it produces a small but distracting gauge pulse each time the horizon crosses a grid zone. 
Horizon drift can be eliminated by introducing a feedback mechanism that adjusts the magnitude of the shift 
to compensate. We do this as follows: If we wish to force the horizon to remain at some radius tq, but we find the 
horizon is actually at radius tah, then we set 

or _N [rAH - rp) 
^ - ' 

where [3^ , N, and A refer to values at the horizon location vah- This feedback mechanism is not necessary for 
sufficiently fine grid spacing. 

While locating an apparent horizon on a numerically generated time slice is a difficult problem in multidimensions 
|p2|-p5|, it is trivial in spherical symmetry. The marginal outer trapped surface equation 

Dis' + s's^K.j - Kl = 0, (39) 

where is the spatial unit normal to the surface, reduces to 

d{r) = TrrlA - 2BrKT = (40) 

for a spherically symmetric system. Here we have used our variables defined in Eqs. (^ and (^). We find the apparent 
horizon by first evaluating i3(r) at each grid point, and then by using 3-point interpolation to locate the outermost 
root that satisfies i?'(r) > 0. When locating the apparent horizon after the Macormack predictor step, we must be 
sure to use the value of r and not f in evaluating the function ■i?(r) at each grid point. 

In principle, one can also use a shift condition at the apparent horizon to move a black hole through the numerical 
grid. This could be accomplished by adjusting the shift so that grid points on one side of the hole fall into the horizon, 
and grid points on the other side emerge from it. A similar idea could in principle be applied to rotating black holes 
or systems with more than one black hole. 
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B. Inner Boundary Condition on Eqs. 



Treating the inner boundary correctly is a primary motivation for using a hyperbolic formulation of the Einstein 
equations. The spherically symmetric CBY equations have only three characteristics at each spacetime point. These 
lie along the ingoing and outgoing null rays, and along the normal vector t. A boundary condition is required only 
at a point that cannot obtain information from one or more of the characteristics passing through it. For example, 
the outer boundary cannot obtain information from the ingoing null characteristics that it intersects, because these 
characteristics originate from outside the computational domain, where we have no data. To update quantities at the 
outer boundary, this information must be provided by a boundary condition. Similarly, the inner boundary ordinarily 
requires a boundary condition because it cannot obtain information from the outgoing null characteristics that it 
intersects. However, if the inner boundary follows an outgoing spacelike or null trajectory, then each point on the 
inner boundary can obtain information from all three characteristics passing through it, so a boundary condition is 
not required. 

This is the essence of an AHBC scheme for treating the inner boundary: by forcing the inner boundary to move 
along with the apparent horizon, we force it along an outgoing non-timelike path. Therefore, there is no need to impose 
an explicit inner boundary condition. Regardless of the mathematical formulation of Einstein's equations being used, 
general relativity tells us that when the inner boundary follows an outgoing spacelike or null path, information with 
physical content cannot penetrate this boundary from the inside, since this information cannot propagate outside the 
light cone. A key advantage of a hyperbolic formulation with only simple (non-spacelike) characteristics is that in 
such a formulation, this statement applies to gauge information as well. 

The way we solve Eqs. without imposing an explicit condition at the inner boundary is by simply ignoring the 
innermost grid point during the interpolation step of our causal differencing scheme. In the case where the inner 
boundary r — r,nin moves with respect to the (f, f) coordinate a distance less than Af during each time step, the 
interpolation becomes an extrapolation at the innermost point. This is always the case in our simulation because of 
the Courant limit: Because the inner boundary, which moves along with the apparent horizon, has a velocity with 
respect to the (i,f) coordinate system of approximately N/A, the inner boundary can never move farther than Af 
during a time step without violating the Courant condition At < {A/N)/^f. One could avoid this extrapolation by 
using an implicit differencing scheme jl^,^ to get around the Courant condition, but such a scheme requires much 
more computer time than explicit schemes, especially in the multidimensional case. 



C. Shift in the Remainder of the Spacetime 



Once a shift criterion at the apparent horizon has been chosen, one must then determine the shift in the remainder of 
the spacetime. At spatial infinity, one presumably would like the shift to approach zero, so that the spacetime metric 
components approach Minkowski values. Or perhaps, in the case of spacetimes with nonzero angular momentum, 
one would like the asymptotic shift to describe a co-rotating frame. However, given the shift at infinity and at the 
apparent horizon, it is not clear how to choose the shift elsewhere. 

One possibility is to use a parameterized analytic function whose parameters are set so that the shift behaves 
appropriately near the apparent horizon and far from the black hole. For example, when evolving a single black hole 
in spherical symmetry we have tried the Gaussian form 

^^^-{r-r,f/w^ ^ (41) 

where C, Tc, and w are chosen so that (3^ is equal to N/A at the apparent horizon, dj3^ /dr is equal to d{N/A)/dr at 
the apparent horizon, and is smaller than some threshold at the outer boundary of the grid. Although this choice 
results in a second order convergent evolution, we find that the grid points tend to compress or stretch where the shift 
gradients are large, and eventually coordinate singularities develop that cause the simulation to terminate. Similarly, 
one can choose 

/3'' = ^(l-tanh((r-r,)/«;)), (42) 

so that the shift is equal to N/A far inside some arbitrary radius r = Tc, and zero far outside r = Vc- In this case, the 
grid points become compressed near r = Tc as one evolves in time, and again the simulation terminates. 

It therefore appears that in addition to a prescription for specifying the shift at the apparent horizon and at infinity, 
one must impose some additional restriction on the shift that ensures that it will not induce any coordinate pathologies 
elsewhere. Such a restriction can be provided by two different elliptic shift conditions that were introduced for the 
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very purpose of minimizing coordinate strain caused by a shift vector. The first is the minimal distortion condition 
lE6f, which can be written in the form 



2D^ 



N K 



(43) 



The minimal distortion shift minimizes the average change of shape of a spatial volume element as it is dragged from 
one time slice to the next. A related choice is the minimal strain condition [E6|, which can be written in the form 



(44) 



The minimal strain shift minimizes the average change in the three metric gij as one evolves from one time slice to 
the next. It differs from the minimal distortion shift in that it takes into account the change in size of spatial volume 
elements as well as their change in shape. 

The downside of these shift conditions is that they require one to solve elliptic equations. This can be costly in 
terms of computer time, especially in three dimensions. It may be possible to use a parameterized analytic function to 
mimic one of these conditions, or it may suffice to use an approximate solution. However, since it appears that these 
conditions give us a useful shift vector, we will adopt them in the spherically symmetric case, where the computational 
burden is not so severe. 



Both the minimal distortion and minimal strain conditions work well in spherical symmetry, as shown in section Vl 



They prevent grid points from becoming locally compressed or stretched to the point where coordinate singularities 
form. Using the variables defined in (0) and (g), we find that the minimal distortion equation (H3) takes the form 



~\ 2 



dr 



and the minimal strain equation 



Br 



d^r \ df d 



dr dr 2 Br 



= N ar iK\ - Kt) + 
3) becomes 



Qf2 



-pr 



-13'' 
dr dr 

TrT \ 

2Br J 



Mr 



- M, 



rT 



'AMtt 



(45) 




dr df 



X rr 



TrT 

2Br 



(46) 



We have written both equations in the f coordinate system and included the factors df/dr and d^f/dr^ because the 
shift must be computed not only after the corrector step of the Macormack scheme when f ~ r, but also after the 
predictor step, when f ^ r. 

Eqs. m3t) and (MS) require boundary conditions at both ends of the numerical grid. We impose the condition 



df d 
dr df 



(n > 1) 



at the outer boundary r = rmax so that the spacetime is asymptotically Minkowski, and we impose 

at the apparent horizon, so that the apparent horizon is stationary with respect to the coordinates. 
Both the minimal distortion and minimal strain equations can be written in the general form 

d 

f3- + 2Q-^p- + P^- = R. 



(47) 



(48) 



(49) 
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We solve this equation using the usual three-point finite difference approximation 

^ - 2/3[ + (3l_,) + - PU) + m = i?.. (50) 

In order to retain second-order accuracy, we must be careful always to impose the boundary condition ([47|) at the 
point r = Tynax , which may be different from the outer boundary of the grid, f = rmax ■ We write Eq. (|47| ) in the 
second-order accurate finite-difference form 

+ (A(A + + 2(1 - X^)P: + A(A - l)/3r_i) = 0, (51) 

where 

= — — • (52) 



1 = 1 



max ■ 



Here fmaa; is the value of r corresponding to the point r = rmaj^, and Eqs. (|5l| ) and ( |5^ ) are evaluated at 
Combining Eqs. (50) and ( pT| ) to eliminate the fictitious grid point i„iaa; + 1 yields a boundary condition for Pl^^^ in 



terms of 

To impose the boundary condition ( [48[ ) at the horizon, we use the finite-difference expression 

A(A + l)/3r+i + 2(1 - X')P: + A(A - l)/3[_i = 2 ^ 

where 



(53) 



A^^^. (54) 
Ar 

Here tah is the f-coordinate location of the apparent horizon, and Eqs. (B3h and (M) are evaluated at the value of i 



such that i — 1 is the innermost grid point that lies outside the horizon. Combining Eqs. ( pO| ) and (53) to eliminate 
the grid point at i -I- 1 yields a boundary condition for in terms of /3[. 

We currently use two methods to solve the linear system of equations that results from Eqs. (50), (^l|), and (|5^). 
We use a standard tridiagonal algorithm [ p7[ when running in serial, or when running in parallel using a small number 
of processors. Although very efficient, this is a serial algorithm, so it begins to dominate the computation time as the 
number of processors increases. For a large number of processors, we instead use a multigrid technique | |2^ , p^ which 
is parallelized with the help of the DAGH system. 

V. OUTER BOUNDARY 

The only variables that require a boundary condition at the outer boundary r = r^ax are the quantities Lrr, Lt, 
Mj-rr, Mj-T, o,Qr, and ttrr, which propagate along the ingoing and outgoing null characteristics. All other variables 

evolve along the characteristics parallel to the vector t, and are never differentiated with respect to r in Eqs. (|35|). 

To apply a boundary condition on these variables, we assume that for large r, the function / for which we need a 
boundary condition behaves like 

Ci ^ (55) 

j,n j,n-\-l ^ ' 

for some integer n and constants Cj . We then impose 

^^(^"/) + 2|:(^"/)-o (56) 



at the outer boundary. This boundary condition truncates the series (55) after the second term, that is, it uses the 
correct values of C\ and C2, but sets Cj to zero for j > 2. We set the value of n for each variable according to its 
analytic falloff rate for the static solution. 
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We call this boundary condition an "extended Robin" method because of its similarity to the familiar Robin 
boundary condition 

-^Kf) = Co, (57) 
which one often imposes on a function that behaves like 

/ = Co + ^ + ---. (58) 

Because Eq. (|5^) involves approximating the large-r behavior of the function /, it introduces errors that are no 
longer second-order in the grid discretization, but fall off rapidly as one increases the radius of the outer boundary. 
These errors result solely from the fact that Eq. ( |56| ) is not exactly satisfied by the initial data. To obtain a second- 
order convergent scheme, we effectively eliminate these approximation errors by increasing rmax until the second-order 
error terms dominate. 

Because we use a cell-centered grid, the outer boundary is located at half a grid spacing beyond the outermost grid 
point: Trnax = ^i„„^+i/2- Therefore, to ensure second-order convergence in Ar, we must impose Eq. ( p6[ ) not at the 
outermost grid point ri^^^, but at the outer boundary Tmax- This is necessary because r^^^^^ depends on the grid 
spacing Ar. 

We impose Eq. ( p6| ) by using the following finite-difference approximations for the first and second derivatives: 

^^^+1/2 = ^1(^(7^-3 - 40K,_2 + 102K,_i - 1121:, + 43i:,+i) 

+ 0{Arf, (59a) 

^>;+i/2 = 7^iY,.2 - 3r,-i - 2iy, + 23r,+i) + OiArf. (59b) 
or ' 24Ar 

Substituting these expressions in Eq. ( p6[ ) with i = imax yields an equation that can be solved for the quantity r"/ 
at the fictitious grid point i — imax + 1 • Using this quantity, we can now evaluate the usual one-sided first derivative 
operator in Eqs. (p7|-p8|) at the outermost grid point i = imax- 

This boundary condition needs to be applied only during the predictor step of the Macormack scheme discussed 
in section III C . This is because the corrector step uses a backwards one-sided first derivative operator that is well- 
defined for the outermost grid point. Therefore, we need not worry about whether we use r or f in Eqs. ( |56| ) and (|5S|), 
because the normal and computational coordinate systems coincide at the beginning of the predictor step. 

The above boundary condition works well in the case of a spherically symmetric, static solution, but it does not 
properly handle wave propagation. This includes not only physical gravitational waves, which are absent in spherical 
symmetry, but also wavelike gauge modes, which are present in our numerical solution. While there exist many 
methods for imposing wavelike outer boundary conditions, we find for this one-dimensional problem that it suffices 
to simply move the outer boundary far from the origin. 

VI. DIAGNOSTIC TESTS 

In this section we present results of rigorous convergence tests performed with our code. 

In all cases below, we evolve a Schwarzschild black hole in spherical symmetry. We begin each simulation with 
time-symmetric initial data corresponding to the v = plane of the Kruskal diagram. The initial three-metric is 
written in isotropic coordinates: 

(^)ds^= + [dr^+r^ide^+sinHdcj)^)], (60) 

where M is the mass of the black hole, so that the initial values of A, B, T^rr, ^rT, Km Kt, Mrrr, and M^t are 
given by 

A = B=(l + -) (61) 



' 1 + — (62) 



r \ 2r J 
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Krr = Kt = Mrrr = MrT = 0. (64) 

With harmonic shcing, one is free to choose the initial value of the lapse function. We choose 

iV=(^l + -j , (65) 
so that the lapse agrees with the canonical Schwarzschild expression at large r. This choice of lapse yields 



1 + (1 + -) (1 + -- — ) (67) 



and the constraint equations (nO|-fill) give us 



AP f MY' f MY f 3M M^\ , , 

M2 / MY^ ( MY^ M 



^- = -1^+2^; \^^Y) \^^-r) (^0) 
Mtt = aor 0. (71) 

We truncate our numerical grid just inside the event horizon at r = A//2, and evolve only the exterior region. 

A. Self-Consistent Convergence 

We have designed our code to be second-order accurate in both space and time, so that if we evolve for a fixed 
time our accumulated truncation error should be ©(Ar)^, where Ar is our grid discretization. Any deviation from 
second-order convergent behavior would indicate a mistake in one of our methods or a coding error. 

To test whether our code is indeed second-order convergent, we evolve the same initial data for the same time t 
using different grid discretizations, and we compare the results. Let X(i; Ar) be the value of some quantity X at time 
t, computed with grid discretization Ar. For second-order convergent behavior, 

X(t;Ar) =X(t)t,uc + (Ar)2x(t) 

error ~ • ■ ■ , (72) 

where X(t)tiue is the true solution, and (Ar)^X(/;)oiior is the leading-order truncation error. Hence 

X{t- ^rl2) - X{t; Ar) = -^{Ar)^X{tY,, + ..., (73a) 

4 [X{t; Ar/4) - X{t; Ar/2)] = -^{ArfX{tYor + ..., (73b) 

16 [X{t; Ar/8) - X{t; Ar/4)] = -^{Ar)^X{tY,, + ..., (73c) 



and so on, so that each of the left-hand sides of Ens. (|7^) are equal to leading order. 

Figures ^ and^ demonstrate that the relation (|7^) holds for our code. These figures show results using five different 
grid resolutions, each differing by a factor of 2. They were produced by evolving to t = 5.625M using our AHBC 
method with a minimal distortion shift. This time corresponds to 240 time steps on the coarsest grid, and 3840 time 
steps on the finest. The outer boundary is placed at r = 6AM. Results from different resolutions are subtracted and 
scaled according to Eq. (|7^), and plotted together. 

If only second-order error terms were present, the four plots in each figure would coincide. This is indeed the case 
except at small r, where the plots differ slightly because of third-order and higher error terms, which are present 
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FIG. 3. Second-order error in the variable ot at time t = 5.625M. The label aT{N) denotes the value of ar computed 
using N grid points per unit r. The apparent horizon is located at approximately r — M/2. All four plots coincide except at 
small values of r, where higher-order error terms become significant. The portions of the graph at small r where ar appears 
to become constant represent grid points that are not included in the evolution because they are inside the apparent horizon. 
The values of ar at these grid points are simply set to the value at the innermost significant point for purposes of the plot. 
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FIG. 4. Second-order error in the variable Lrr at time t — 5.625A/, for the same case as shown in Figure |^. 

but not explicitly given in Eqs. (fz^). These error terms vanish faster than (Ar)^ as Ar — > 0, as indicated by the 
convergence of the plots toward each other as the grid resolution is increased. 
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FIG. 5. Second-order error in the variable r (3^ at time t = 5.625Af, for the same case as shown in Figure g. 

Figure |^ shows the second-order error in the quantity r'^(3^ , for the same case shown in Figures I and |. This 
demonstrates that our solution of the minimal distortion equation is second-order convergent, and that it satisfies the 
outer boundary condition to 0{Ar)^. Higher-order error terms are also present in this figure. They converge to zero 
as the grid resolution is increased. 

B. Convergence of Constraint Equations 



Even if a numerical code is second-order convergent in the self-consistent manner described above, it still may not 
converge to the analytic solution. To test whether this is the case, we evaluate the left-hand sides of the constraint 
equations (pO|-pT|), which are not explicitly solved in the code except on the initial time slice. If we are indeed solving 
Einstein's equations, these quantities should vanish like (Ar)^. 

Figure ^ shows left-hand side of Eq. (lOt) for the same case shown in Figures ^-||, and for five different grid 
resolutions. For each finer grid discretization, this quantity is multiplied by a factor of 4, so that for strict second- 
order convergence, the five plots should coincide. Because they do coincide except for small r where higher-order error 
terms are significant, we see that the momentum constraint is satisfied to 0{Ar)^. The same is true for the other 
constraints. 

There is an additional test we can use to determine whether we are solving the correct equations: For a static black 
hole, the total mass inside each radius r, which is a well-defined locally measurable quantity in spherical symmetry, 
should be conserved. The degree to which mass conservation is violated provides a good indicator for the overall 
accuracy of the code. An invariant expression for the mass inside a spherical surface of symmetry is ||29(| 



M 



A 



1 



16nA 



where A is the area of the surface. In terms of our variables, this expression reduces to 



M = 



Br 



1 - 



2A 



(74) 



(75) 



Figure ^ shows the quantity 
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FIG. 6. Momentum constraint at time t = 5.625M, for the same case as shown in Figure ^. In the plot labels, J^{N) denotes 
the value of the left-hand side of Eq. (10b) computed using N grid points per unit r. 
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FIG. 7. The quantity AM at time t — 5.625Af, for the same case as shown in Figure |^. 
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at time t — 5.625M as calculated from Eq. (|75|), for five grid discretizations. This quantity is multiplied by factors of 
4 for each finer grid discretization. The mass of the system is conserved to 0{Ar)'^. 
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VII. DISCUSSION 



While the CBY formulation of Einstein's equations at first glance looks more complicated than the usual ADM 
formulation, in most respects it is actually much simpler, particularly from a numerical point of view. Each equation in 
the CBY system is either a wave equation with an advective term, or a simple advective equation. In the coordinate 
system that is used for causal differencing, the equations are even simpler: They are either wave equations or first- 
order ordinary differential equations in time. Admittedly, the right-hand sides of these equations are complicated 
and nonlinear, but because these right-hand sides contain no derivatives, they do not make numerical solution of the 
equations any more difficult. The large number of evolution equations in the CBY formalism is also not a serious 
drawback because the equations are all of the same form, so they can all be solved by the same method. 

Unlike a previous work |^ in which we forced the horizon to sit at a particular grid point, here we allow the 
horizon to lie at some arbitrary location on the grid. This is a closer approximation to what we expect in the three- 
dimensional case, where one would most likely use Cartesian coordinates to describe a spherical hole. Likewise, in 
the same work we imposed explicit boundary conditions on the apparent horizon in order to solve elliptic constraint 
equations outside the hole. Here we do not impose explicit boundary conditions (except for the shift via Eq. (|4^)) 
because in the three-dimensional case it is not only difficult to impose such conditions on a nonspherical boundary, 
but also the number of boundary conditions needed for solving all of the constraints in three dimensions exceeds the 
number of boundary conditions available at the horizon. 

A key milestone in three-dimensional black hole simulations is the ability to stably move a hole through the numerical 
grid. This is arguably a necessary precursor to simulating binary orbits. Techniques for moving holes through the 
grid could be tested in spherical symmetry using our code, and we plan such tests in the future. To preserve spherical 
symmetry, the motion in this case would be the expansion or contraction (in coordinate space) of the hole rather than 
the translation of the hole's center, but since the location of the center is irrelevant for AHBC methods, many of the 
same methods should apply to both cases. 

The code as described here tends to lose accuracy at times greater than 10 or 20 M, where M is the mass of the 
black hole. Because this behavior is relatively independent of numerical details such as grid resolution, we believe 
that it may be due to either gauge modes or unphysical rapidly-growing solutions of the evolution equations that do 
not satisfy the constraints. One way of solving this problem is to enforce several of the constraint equations after 
each time step. With this modification, we can integrate past t = lOOOM. However, the details of constraint- violating 
solutions, gau ge modes, and methods for dealing with them are beyond the scope of this paper and will be dealt with 
elsewhere IitI. 
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